function [parameters] = dGemricT1fit(serie,masknr, sigma, basedir)

% Estimates parameters and errorbars for T1 recovery in an IR sequence 
% model = A - B*exp(-TI/T1) with dGEMRIC folder structure

% The estimated parameters are saved in a parameters.mat:
% parameters[parameter,serie,prot,masknr,mask] with parameter code:
% 1: meanT1 2:medianT1 3:stdT1 4:weightedmeanT1
% 5:weightedSqMeanT1 6:medianCR

% Output are dicom images of the Amap, Bmap, R1map, R1errormap, T1map and T1errormap, saved
% in *seriedir*/postprocessing/maps 

% IMPORTANT: The values of the R2 and R2 errormap are multiplied with 10000 to get into sensible gray values.
% The errors in the errormap are the lower bounds on the standard deviations calculated with
% the Cramer-Rao lower bounds.

% Henk Smit, h.smit@erasmusmc.nl, 10-2011

cd(basedir)
if exist('locations.mat')
    load('locations.mat')
else
    error( 'dGemricT1fit:nolocations', 'No locations.mat found in basedir. First draw rois for this serie!');  
end

nrimages = 5;%nr of inversion times
protocols = 2;%size(locations(serie).loc,1);
CRLBsigma = sigma;
masknames={'_fem_post','_fem_wb' '_plat_wb'};
protocolnames={'IR_2100_','non_ARC_','Qmap_ARC_','High_Res_'}; %REVERSED NON ARC AND QMAP ARC!!

cd(locations(serie).loc(1).dir)
cd ..
volunteerdir = pwd;
d = dir;
cd ..
postprocdir = fullfile(pwd,'Postprocessing');
qmapdir = locations(serie).loc(2).dir;


% read in data in structure: fitdata(protocol,image)
% and masks in structure: masks(protocol,slice,mask)
disp(['Reading in dicom images']);
for(prot=1:protocols)
    directory = locations(serie).loc(prot,masknr).dir;
    file=locations(serie).loc(prot,masknr).file;
    slidenr(prot)=str2num(file(1:4));
    if prot==1 %Read In SPGRs
            k=1; %image counter
            bone=1; %femur/tibia
            for j=3:size(d,1)
                x=d(j).name;
                if size(char(regexp(x,'TIR_2100', 'match')),2)>=1 %read in the fixed longest IR image
                    k=1;
                    info = dicominfo([volunteerdir, filesep,d(j).name,filesep,file]);   
                    fitdata(prot,k,1).image = double(dicomread(info));
                    fitdata(prot,k,1).TE=info.EchoTime;
                    fitdata(prot,k,1).TR=info.RepetitionTime+info.InversionTime;
                    fitdata(prot,k,1).Angle=info.FlipAngle;
                    fitdata(prot,k,1).info=info;
                    fitdata(prot,k,1).TI=info.InversionTime;
                    fitdata(prot,k,2).image = double(dicomread(info));
                    fitdata(prot,k,2).TE=info.EchoTime;
                    fitdata(prot,k,2).TR=info.RepetitionTime+info.InversionTime;
                    fitdata(prot,k,2).Angle=info.FlipAngle;
                    fitdata(prot,k,2).info=info;
                    fitdata(prot,k,2).TI=info.InversionTime;
                end   
            end
                
            cd(fullfile(postprocdir,'Registrations','femur'));
            dregis=dir;
            for j=3:size(dregis,1) % Read in the registered SPGR images
                k=k+1;
                bone=1;
                info = dicominfo([dregis(j).name,filesep,file]);    
                fitdata(prot,k,bone).image = double(dicomread(info));
                fitdata(prot,k,bone).TE=info.EchoTime;
                fitdata(prot,k,bone).TR=info.RepetitionTime+info.InversionTime;
                fitdata(prot,k,bone).Angle=info.FlipAngle;
                fitdata(prot,k,bone).info=info;
                fitdata(prot,k,bone).TI=info.InversionTime;
            end       
            
            cd(fullfile(postprocdir,'Registrations','tibia'));
            dregis=dir;
            k=1;
            for j=3:size(dregis,1) % Read in the registered SPGR images
                k=k+1;
                bone=2;
                info = dicominfo([dregis(j).name,filesep,file]);    
                fitdata(prot,k,bone).image = double(dicomread(info));
                fitdata(prot,k,bone).TE=info.EchoTime;
                fitdata(prot,k,bone).TR=info.RepetitionTime+info.InversionTime;
                fitdata(prot,k,bone).Angle=info.FlipAngle;
                fitdata(prot,k,bone).info=info;
                fitdata(prot,k,bone).TI=info.InversionTime;
            end

    else     % Read in QMAPS
        for i=1:nrimages
%             qmapregdir = fullfile(volunteerdir, 'Postprocessing','RegistrationsQMAP',directory(end-8:end-1));
            filenr=slidenr(prot)+4*(i-1);
            qbone = 1;
            if filenr>9
                info = dicominfo(fullfile(qmapdir,['00',num2str(filenr),'.dcm']));
            else
                info = dicominfo(fullfile(qmapdir,['000',num2str(filenr),'.dcm']));
            end
            fitdata(prot,i,qbone).image = double(dicomread(info));
            fitdata(prot,i,qbone).image = double(dicomread(info));
%             fitdata(prot,i,qbone).TE=info.EchoTime;
            fitdata(prot,i,qbone).TR=info.InversionTime+info.RepetitionTime;
            fitdata(prot,i,qbone).Angle=info.FlipAngle;
            fitdata(prot,i,qbone).info=info;
            fitdata(prot,i,qbone).TI=info.InversionTime;
        end
    end

    % Read in masks: masks(prot,masknr,mask) mask:1=femp,2=femwb,3=plat

        
        cd(fullfile(postprocdir,'Mapping_Masks'));
        dmask=dir;
    for j=3:size(dmask,1)
        m=0;
        x=dmask(j).name;
        for mask=3:size(dmask,1) 
            if size(char(regexp(dmask(mask).name,[protocolnames{prot} file(1:4)], 'match')),2)>=1
                m=m+1;
                load(dmask(mask).name,'M')
                tempmasks(m).mask=M;
            end
        end
        for r=1:3
            masks(prot,masknr,r).mask=tempmasks(r).mask;
        end
        clear tempmasks
    end
end

%conversion factor from standard to rayleigh distribution
% CRLBsigma=1.527*maxsd

clear i j k m r prot mask;
for prot=1:protocols
        for mask=1:size(masks,3)

            disp(['Estimating map on mask ' num2str(mask) '/' num2str(size(masks,3)) ' on slice: ' num2str(masknr) '/' num2str(size(masks,2)) ' on protocol ' num2str(prot) '/' num2str(protocols) ' on serie ' num2str(serie)]);

            for img=1:nrimages
                if mask<3 || prot==2
                    maskedfitdata(prot,masknr,mask,img).image=fitdata(prot,img,1).image.*masks(prot,masknr,mask).mask; %Femur
                elseif mask==3 && prot==1
                    maskedfitdata(prot,masknr,mask,img).image=fitdata(prot,img,2).image.*masks(prot,masknr,mask).mask; %Tibia
                end
            end

            if prot==1
                for img=1:nrimages
                    TI(img,1)=double(fitdata(prot,img,1).TI);
                    TR(img,1)=double(fitdata(prot,img,1).TR);
                end
            else
                for img=1:nrimages
                    TI(img,1)=double(fitdata(prot,img,1).TI);
                    TR(img,1)=double(fitdata(prot,img,1).TR);
                end
                
            end

            yd=zeros(size(TI,1),size(TI,2));

            [row,column]=find(maskedfitdata(prot,masknr,mask,img).image);
            nrvoxels=size(row,1);


            dims=size(maskedfitdata(prot,masknr,mask,img).image);
            T1=zeros(dims(1),dims(2));
            R1=T1;
            A=T1;
            B=T1;
            STDCRT1=T1;
            STDCRR1=T1;
            STDCRB=T1;
            STDCRA=T1;
            yd=zeros(nrimages,1);

        opt = optimset('fminunc');
        opt = optimset(opt,'Diagnostics','off','DerivativeCheck','off','LargeScale','off','gradObj','on','Display','off','MaxIter',200,'Hessian','off','TolFun',1e-12,'Tolx',1e-10,'MaxFunEvals',200);
        s = fitoptions('Method','NonLinearLeastSquares','Lower',[0,0,0],'Upper',[30000,30000,100],'Startpoint',[5000, 1.3, 0.001],'Maxiter',20,'Display','off','TolFun',10^-5,'TolX',10^-5);
        % f = fittype('a-(b*exp(-x*c))','options',s);
        fiteq = (['a*(1-b*exp(-x*c)+ exp(-(' num2str(TR(1)-TI(1)) '+x)*c))' ]);
        f = fittype(fiteq,'options',s);
        fitrange = 2:size(TI);

        for i=1:nrvoxels
    %         disp(['Calculating pixel: ' num2str(i) '/' num2str(nrvoxels(1))]);
                    if maskedfitdata(prot,masknr,mask,img).image(row(i),column(i)) > 0;
                        for k=1:nrimages
                            yd(k,1)=double(maskedfitdata(prot,masknr,mask,k).image(row(i),column(i)));
                        end

                        [LS,ML] = T1IRAbsComputePar(yd,TI,0,CRLBsigma,opt,s,f,TR,fitrange);
                        [CR] = T1IRAbsCramerRao(ML,TI,CRLBsigma,TR);

                        A(row(i),column(i))=ML(1,1);
                        B(row(i),column(i)) = ML(2,1);
                        T1(row(i),column(i)) = 1/ML(3,1);
                        R1(row(i),column(i)) = ML(3,1);

                        if CR(1,1)>0 && ~isnan(CR(1,1))
                            STDCRA(row(i),column(i)) = sqrt(CR(1,1));
                        else
                            STDCRA(row(i),column(i)) = 10000;
                        end

                        if(CR(2,2)>0 && ~isnan(CR(2,2)))
                            STDCRB(row(i),column(i)) = sqrt(CR(2,2));
                        else
                            STDCRB(row(i),column(i)) = 10000;
                        end

                        if(CR(3,3)>0 && ~isnan(CR(3,3)))
                            STDCRR1(row(i),column(i)) = sqrt(CR(3,3));
                            STDCRT1(row(i),column(i)) = sqrt(CR(3,3))/(ML(3,1)^2); %error propagation dT/T = dR/R
                        else
                            STDCRR1(row(i),column(i)) = 5;
                            STDCRT1(row(i),column(i)) = 10000;
                        end

                    else
                        A(row(i),column(i))= 0;
                        B(row(i),column(i)) = 0;
                        T1(row(i),column(i)) = 0;
                        R1(row(i),column(i)) = 0;

                        STDCRR1(row(i),column(i)) = 0;
                        STDCRT1(row(i),column(i)) = 0;
                    end  

            end

            %write result in dcm files. R1 values*10000 to get in a sensible range
            if slidenr(prot)>9
                    filenr = ['00',num2str(slidenr(prot))];
            else
                    filenr = ['000',num2str(slidenr(prot))];
                end
            outputdir = fullfile(postprocdir,'Maps');
            file = slidenr(prot);
            outputdirfile = [outputdir,filesep,filenr,masknames{mask},protocolnames{prot}];
            if ~exist(outputdir, 'dir')
              disp( ['making: ' outputdir]);     
              [success, message, messageid] = mkdir(outputdir);
              if ~success    
                 error( 'ratheartT2fit:invalidoutputdir', message );     
              end 
            end
                dicomwrite(int16(10000*R1),[outputdirfile,'_R1map','.dcm'],fitdata(prot,1).info);
                dicomwrite(int16(T1),[outputdirfile,'_T1map','.dcm'],fitdata(prot,1).info);
                dicomwrite(int16(A),[outputdirfile,'_Amap','.dcm'],fitdata(prot,1).info);
                dicomwrite(int16(B),[outputdirfile,'_Bmap','.dcm'],fitdata(prot,1).info);
                dicomwrite(int16(10000*STDCRR1),[outputdirfile,'_R1errormap','.dcm'],fitdata(prot,1).info);
                dicomwrite(int16(STDCRT1),[outputdirfile,'_T1errormap','.dcm'],fitdata(prot,1).info);

                if exist([basedir filesep 'Estimated_parameters.mat'])==2
                    load([basedir filesep 'Estimated_parameters.mat'])
                end            

                %1: meanT1 2:medianT1 3:stdT1 4:weightedmeanT1
                %5:weightedSqMeanT1 6:medianCR

                parameters(1,serie,prot,masknr,mask) = mean(nonzeros(T1));
                parameters(2,serie,prot,masknr,mask) = median(nonzeros(T1));
                parameters(3,serie,prot,masknr,mask) = std(nonzeros(T1));
                parameters(4,serie,prot,masknr,mask) = (nonzeros(T1)'*(1./nonzeros(STDCRT1)) ) / sum(1./nonzeros(STDCRT1)); %T1/error
                parameters(5,serie,prot,masknr,mask) = (nonzeros(T1)' * (1./(nonzeros(STDCRT1).^2)) ) / sum(1./(nonzeros(STDCRT1)).^2); %T1/error^2
                parameters(6,serie,prot,masknr,mask) = median(nonzeros(STDCRT1));

                save([basedir, filesep,'Estimated_parameters.mat'],'parameters')
    %             imshow(T1,[]);
        end
%     end
end
    
